module util_cmfd_eddington
    implicit none
!    real(8), pointer, dimension(:,:) :: hDif, hTot, hRmv, hNuF
!    real(8), pointer, dimension(:,:,:) :: hmSct
    real(8), private, pointer, dimension(:,:) :: mm2, psa
    real(8), pointer, dimension(:,:) :: hE
    integer, private, parameter :: fdhat = 777, fdtilde = 877
    integer, private, parameter :: fj=91, fdif=92, fphiac=93, fcmfd_EDT=1094, fMTCMFD=1095, fMTCMFDOFF=1096, fRL1CMFD=1097, fRL2CMFD=1098
    integer, private, parameter :: fcj=1091, fcp=1093
    real(8), private, pointer, dimension(:, :) :: phi
    real(8), private, pointer, dimension(:) :: psi, psip
    real(8) :: l2ref, sumref
contains

subroutine prepare_cmfd_EDT
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "accm.FI"
    integer :: g, i

    allocate(hE(totm_cmfd, ng))
!    allocate(hTot(totm_cmfd, ng))
!    allocate(hRmv(totm_cmfd, ng))
!    allocate(hNuF(totm_cmfd, ng))
!    allocate(hmSct(ng,ng, totm_cmfd))
!    allocate(phiftemp(totmesh, ng))
    allocate(phi(totm_cmfd, ng))
    allocate(psi(totm_cmfd))
    allocate(psip(totm_cmfd))

!    open(fMTCMFD, file='MT_CMFD_ED.txt', status='unknown')
!    open(fMTCMFDOFF, file='MT_CMFDOFF_ED.txt', status='unknown')
!    open(fRL1CMFD, file='MT_R1CMFD_ED.txt', status='unknown')
!    open(fRL2CMFD, file='MT_R2CMFD_ED.txt', status='unknown')
    phi = 1.
    
    if (PDFREF) then
        sumref = sum(refc)
        l2ref = 0
        do g=1, ng
            do i=1, totm_cmfd
                l2ref = l2ref + refc(i, g)**2
            end do
        end do
        l2ref = sqrt(l2ref)
    end if 
    
!    open(fcmfd_EDT, file='CMFD_EDT.txt', status='unknown')

!    if (NMVWND .eq. 1) then
        mm2 => mmt2nd
        psa => phi_sa
!    else
!        mm2 => mmt2ndAC
!        psa => phisaAC
!    endif

end subroutine

subroutine evaluateDtilde
    implicit none
    
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "hmgz.FI"
    
    integer i, j, k, g
    real(8)    betaL, betaR, compL, compR
    
    do g=1, ng
        do i=1, totm_cmfd
            hE(i,g) = 1._8/hTot(i,g)*hphi2(i,g)/hphi(i,g)
        end do
    end do
    Dtilde_c = 0.
    do g= 1, ng    
        do i=1, totm_cmfd
            ! for left boundary
            if (i.eq.1) then
                compR = i
                betaR = hE(compR,g)/h_cmfd(i)
                Dtilde_c(1, g) = 2*albedo(1)*betaR/(albedo(1)+2*betaR)
!                Dtilde_c(1, g) = 0.
                compL = i
                compR = i+1
                betaL = hE(compL, g)/h_cmfd(i)
                betaR = hE(compR, g)/h_cmfd(i+1)
                Dtilde_c(2, g) = 2*betaL*betaR/(betaL+betaR)
            elseif(i.eq.totm_cmfd) then
                compL = i-1
                compR = i
                betaL = hE(compL, g)/h_cmfd(i-1)
                betaR = hE(compR, g)/h_cmfd(i)
                Dtilde_c(2*totm_cmfd-1, g) = 2*betaL*betaR/(betaL+betaR)
                compL = i
                betaL = hE(compL, g)/h_cmfd(i)
                Dtilde_c(2*totm_cmfd, g) = 2*albedo(2)*betaL/(albedo(2)+2*betaL)
!                Dtilde_c(2*totm_cmfd, g) = 0.
            else
                compL = i-1
                compR = i
                betaL = hE(compL, g)/h_cmfd(i-1)
                betaR = hE(compR, g)/h_cmfd(i)
                Dtilde_c(2*i-1, g) = 2*betaL*betaR/(betaL+betaR)
                compL = i
                compR = i+1
                betaL = hE(compL, g)/h_cmfd(i-1)
                betaR = hE(compR, g)/h_cmfd(i)
                Dtilde_c(2*i , g) = 2*betaL*betaR/(betaL+betaR)
            end if
        end do
    end do
end subroutine

!subroutine condense(add, wnd, wrt)
!    use define
!    use mvwnd
!    implicit none
!    include "pspec.FI"
!    include "arrays.FI"
!    include "accm.FI"
!    include "xsec.FI"
!    integer :: g, gp, i, i_f, i_c
!    real(8) :: rtot, rrmv, rsct(ng), volflux, rvolflux
!    logical,intent(in) :: add, wnd, wrt
!    
!    ! flux condense
!    phi_c = 0.
!    phi_cp = 0.
!    hTot=0;hRmv=0;hmSct=0;hNuF=0
!    
!    ! cycle flux are used to obtain initial guess
!    do g=1, ng
!        do i=1, totmesh
!            i_c = fi2co(i)
!            phi_cp(i_c, g) = phi_cp(i_c, g) + phi_f(i,g)*h(i)
!        end do
!    end do
!    do g=1, ng
!        do i=1, totm_cmfd
!            phi_cp(i, g) = phi_cp(i, g)/h_cmfd(i)
!        end do
!    end do
!    ! surface current condense    
!    do g=1, ng
!        i_f = 1
!        do i=1, totmesh*2
!            if (jcvt2cmfd(i)) then
!                jc(i_f, g) = jf(i, g)
!                i_f = i_f + 1
!            end if
!        end do
!    end do
!    if (wrt) then
!        write(fcp, '(2000(e11.5,1x))') ((phi_cp(i,g),i=1,totm_cmfd), g=1, ng)
!        write(fcj, '(2000(e11.5,1x))') ((jc(i,g),i=1,totm_cmfd*2), g=1, ng)
!    endif 
!    
!    if (add) then
!        call add_mvwnd(phi_f, jc)
!    endif
!
!    if (wnd) then
!        call get_mvwnd(phiftemp, jc)
!    else
!       ! surface current condense    
!        do g=1, ng
!            i_f = 1
!            do i=1, totmesh*2
!                if (jcvt2cmfd(i)) then
!                    jc(i_f, g) = jfAC(i, g)
!                    i_f = i_f + 1
!                end if
!            end do
!        end do
!        phiftemp = phifAC
!    end if 
!    
!     
!    do g=1, ng
!!        rtot = 0;rrmv=0;rsct=0;
!        do i=1, totmesh
!            i_c = fi2co(i)
!            
!            ! accumulated values are used for constructing the system
!            volflux = phiftemp(i,g)*h(i)
!            phi_c(i_c, g) = phi_c(i_c, g) + volflux
!            hTot(i_c,g) = hTot(i_c,g) + xTot(complist(i),g)*volflux
!            hRmv(i_c,g) = hRmv(i_c,g) + xRmv(complist(i),g)*volflux
!            hNuF(i_c,g) = hNuF(i_c,g) + xNuf(complist(i),g)*volflux
!            do gp=1, ng
!                hmSct(g,gp,i_c) = hmSct(g,gp,i_c) + XsctM(g,gp,complist(i))*volflux
!            enddo            
!        end do
!    end do
!    do g=1, ng
!        do i=1, totm_cmfd
!            if (phi_c(i, g)>0.) then
!                rvolflux = 1./phi_c(i, g)
!                hTot(i,g)=hTot(i,g)*rvolflux
!                hRmv(i,g)=hRmv(i,g)*rvolflux
!                hNuF(i,g)=hNuF(i,g)*rvolflux
!!#ifdef _LANISO
!!            hDif(i,g)=1./(3*(hTot(i,g)-0.1_8))            
!!#else
!!            hDif(i,g)=1./(3*hTot(i,g))
!!#endif
!                hDif(i,g)=1./(3*hTot(i,g))
!                do gp=1, ng
!                    hmSct(g,gp,i) = hmSct(g,gp,i)*rvolflux
!                end do
!                phi_c(i, g) = phi_c(i, g)/h_cmfd(i)
!            endif
!        end do
!    end do
!    
!
!    ! generate fission source
!    psi_cp = 0.
!    do g=1, ng
!        do i=1, totm_cmfd
!            psi_cp(i) = psi_cp(i) + hNuf(i, g)*phi_cp(i, g)
!        end do
!    end do        
!end subroutine


subroutine evaluateDhat
    implicit none
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "hmgz.FI"
    integer i, g

    Dhat = 0.
    do g= 1, ng    
        do i=1, totm_cmfd
            ! for left boundary
            if (hphi(i, g) .eq. 0._8) cycle
            if (i.eq.1) then 
                if (hj(i*2-1, g)==0. .and. hphi(i, g)==0) then
                    Dhat(i*2-1, g) = Dtilde_c(i*2-1, g)
                else
                    Dhat(i*2-1, g) = -hj(i*2-1, g)/hphi(i, g)
                end if
                Dhat(i*2, g) = -(hj(i*2, g)+Dtilde_c(i*2, g)*(hphi(i+1, g)-hphi(i, g)))/(hphi(i+1, g)+hphi(i, g));
            elseif(i.eq.totm_cmfd) then
                Dhat(i*2-1, g) = -(hj(i*2-1, g)+Dtilde_c(i*2-1, g)*(hphi(i-1, g)-hphi(i, g)))/(hphi(i-1, g)+hphi(i, g));
                if (hj(i*2, g)==0. .and. hphi(i, g)==0) then
                    Dhat(i*2, g) = Dtilde_c(i*2, g)
                else
                   Dhat(i*2, g) = -hj(i*2, g)/hphi(i, g)
                end if
            else
                Dhat(i*2-1, g) = -(hj(i*2-1, g)+Dtilde_c(i*2-1, g)*(hphi(i-1, g)-hphi(i, g)))/(hphi(i-1, g)+hphi(i, g));
                Dhat(i*2, g) = -(hj(i*2, g)+Dtilde_c(i*2, g)*(hphi(i+1, g)-hphi(i, g)))/(hphi(i+1, g)+hphi(i, g));
            end if
        end do
    end do
    cntacc = cntacc+1
end subroutine

subroutine make_matrix
    use LU_tdiag
    implicit none    
    
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "hmgz.FI"
    integer i, j, k, l, g, left, right
    real(8) ::offsum
    real(8), dimension(totm_cmfd*2) :: Dhatavg
    
    do g=1, ng
        !Dhatavg = Dhatacc(:,g)/cntacc
        Dhatavg = Dhat(:,g)
        mt_cmfd(1, 1, g) = 0.
        mt_cmfd(3, 1, g) = -(Dtilde_c(2, g)+Dhatavg(2))
        mt_cmfd(2, 1, g) = -Dhatavg(1)-Dhatavg(2)+Dtilde_c(2, g)+hRmv(1, g)*h_cmfd(1)
        do i=2, totm_cmfd-1
            left = i*2-1
            right = i*2
            mt_cmfd(1, i, g) = -(Dtilde_c(left, g)+Dhatavg(left))
            mt_cmfd(3, i, g) = -(Dtilde_c(right, g)+Dhatavg(right))
            mt_cmfd(2, i, g) = -Dhatavg(left) - Dhatavg(right) + Dtilde_c(left, g) + Dtilde_c(right,g) + hRmv(i, g)*h_cmfd(i)
        end do
        mt_cmfd(1, totm_cmfd, g) = -(Dtilde_c(totm_cmfd*2-1, g)+Dhatavg(totm_cmfd*2-1))
        mt_cmfd(3, totm_cmfd, g) = 0.
        mt_cmfd(2, totm_cmfd, g) = -Dhatavg(totm_cmfd*2-1)-Dhatavg(totm_cmfd*2)+Dtilde_c(totm_cmfd*2-1, g) + hRmv(i, g)*h_cmfd(totm_cmfd)
    end do
    
    
    do g=1, ng
        call LU_factorize(totm_cmfd, mt_cmfd(:,:,g))
    end do
end subroutine


subroutine make_matrix_withoutcurrent
    use LU_tdiag
    implicit none    
    
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "hmgz.FI"
    integer i, j, k, l, g, left, right
    real(8) ::offsum
    real(8), dimension(totm_cmfd*2) :: Dhatavg
    
    do g=1, ng
        do i=1, totm_cmfd
            left = (i-1)*2+1
            right = i*2
            mt_cmfd(1, i, g) = -Dtilde_c(left, g)
            mt_cmfd(3, i, g) = -Dtilde_c(right,g)
            mt_cmfd(2, i, g) =  Dtilde_c(left, g) + Dtilde_c(right,g) + hRmv(i, g) * h_cmfd(i)
        end do
        mt_cmfd(1, 1, g) = 0.
        mt_cmfd(3, totm_cmfd, g) = 0.
    end do
    do g=1, ng
        call LU_factorize(totm_cmfd, mt_cmfd(:,:,g))
    end do
end subroutine

subroutine makeFS
    implicit none
    include "pspec.FI"
    include "xsec.FI"
    include "hmgz.FI"
    integer :: i, g
    psi =0
    do g=1, ng
        do i=1, totm_cmfd
            psi(i) = psi(i) + hNuF(i,g)*phi(i,g)*h_cmfd(i);
        enddo
    enddo
end subroutine

subroutine makeQ(g, k, q)
    implicit none
    include "pspec.FI"
    include "xsec.FI"
    include "hmgz.FI"
    integer,intent(in) :: g
    real(8),intent(in) :: k
    real(8), dimension(:), intent(out) :: q
    integer :: i,  gp
    q =0
    do i=1, totm_cmfd
        q(i) = Xchi(1,g)*psi(i)/k
        do gp = 1, ng
            if (gp.eq.g)  then 
                continue
            else
                q(i) = q(i) + hmSct(g,gp,i)*phi(i,gp)*h_cmfd(i)
            end if 
        end do
    end do
end subroutine

subroutine prolongation
    implicit none
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    real(8), dimension(totm_cmfd) :: ratio
    real(8) :: rtnorm
    integer :: g, i, i_f, i_c
    
!    rtnorm = sum(psi_f)/sum(psi)
!    do g=1, ng
!        do i=1, totm_cmfd
!            phi(i,g) = phi(i,g) * rtnorm
!            if(phi_c(i,g)> 0) then
!                ratio(i) = phi(i,g)/phi_c(i,g)
!            else
!                ratio(i) = 1.
!            end if
!            
!            if (ratio(i)>2) then
!                ratio(i) = 2.
!            end if
!        end do
!        !ratio = phi_c(:, g)/phi_cp(:, g) * sum(phi_cp(:, g))/sum(phi_c(:,g))
!       
!        do i=1, totmesh
!            i_c = fi2co(i)
!            phi_f(i, g) = phi_f(i, g) * ratio(i_c)
!!            adjp(i, g) = ratio(i_c)
!            
!        end do
!    end do
    do i=1, totmesh
        i_c = fi2co(i)
        if (psi_c(i_c)>0._8) then
            adjp(i) = psi(i_c)/psi_c(i_c)
        else 
            adjp(i) = 1.
        end if 
    end do
    
!    adjp = psi_c/psi_cp

!    ! update fission source
!    psi_f = 0.
!    do g=1, ng
!        do i=1, totmesh
!            psi_f(i) = psi_f(i) + Xnuf(complist(i), g) * phi_f(i, g) * h(i);
!        enddo
!    enddo
    
end subroutine       

subroutine err_coarse(phic, err1, err2, err3) 
    include "pspec.FI"
    include "xsec.FI"
    include "arrays.FI"
    include "hmgz.FI"
    real(8), intent(out) :: err1, err2, err3
    real(8), intent(in), pointer, dimension(:, :) :: phic
    integer :: g, i
    real(8) :: rn ! ratio of normalization
    real(8) :: dif
    if (PDFREF .eq. .false.) return

    err1 = 0;err2 = 0;err3 = 0
    rn = sum(phic)/sumref
    do g=1, ng
        do i=1, totm_cmfd
            dif = abs(refc(i,g)*rn-phic(i,g))
            dif = dif /(refc(i,g)*rn)
            err1 = err1 + (dif)**2
            if (dif > err2) err2 = dif
            if (hNuF(i, g)>0 .and. dif>err3) err3 = dif
        end do
    end do
end subroutine


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!         W i e l a n d t         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine make_matrix_full
    use LU_general
    implicit none    
    
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "hmgz.FI"

    integer i, j, k, l, g, left, right, gp
    real(8) ::offsum
    real(8), dimension(totm_cmfd*2) :: Dhatavg
    
    mtfull = 0.
    do g=1, ng
        !Dhatavg = Dhatacc(:,g)/cntacc
        Dhatavg = Dhat(:,g)
        i=1
        l=(g-1)*totm_cmfd + i
        !mtfull(l, l) = -Dhatavg(i)-Dhatavg(i+1)+Dtilde_c(i+1, g)+hRmv(i, g)*h_cmfd(i)
        mtfull(l, l) = -Dhatavg(i)+Dtilde_c(i+1, g)+hRmv(i, g)*h_cmfd(i)
        do gp=1, ng
            if (gp .ne. g) mtfull(i+(gp-1)*totm_cmfd, l) = -hmSct(g, gp, i)*h_cmfd(i)
        end do
        !mtfull(l+1, l) = -(Dtilde_c(i+1, g)+Dhatavg(i+1))
        mtfull(l+1, l) = -(Dtilde_c(i+1, g))
        do i=2, totm_cmfd-1
            l=(g-1)*totm_cmfd + i
            left = i*2-1
            right = i*2
!            mtfull(l-1, l) = -(Dtilde_c(left, g)+Dhatavg(left))
            mtfull(l-1, l) = -(Dtilde_c(left, g))
!            mtfull(l, l) = -Dhatavg(left) - Dhatavg(right) + Dtilde_c(left, g) + Dtilde_c(right,g) + hRmv(i, g)*h_cmfd(i)
            mtfull(l, l) = Dtilde_c(left, g) + Dtilde_c(right,g) + hRmv(i, g)*h_cmfd(i)
            do gp=1, ng
                if (gp .ne. g) mtfull(i+(gp-1)*totm_cmfd, l) = -hmSct(g, gp, i)*h_cmfd(i)
            end do
!            mtfull(l+1, l) = -(Dtilde_c(right, g)+Dhatavg(right))
            mtfull(l+1, l) = -(Dtilde_c(right, g))
        end do
        l=(g-1)*totm_cmfd + i
!        mtfull(l-1, l) = -(Dtilde_c(i*2-1, g)+Dhatavg(i*2-1))
        mtfull(l-1, l) = -(Dtilde_c(i*2-1, g))
!        mtfull(l, l) = -Dhatavg(i*2-1)-Dhatavg(i*2)+Dtilde_c(i*2-1, g) + hRmv(i, g)*h_cmfd(i)
        mtfull(l, l) = -Dhatavg(i*2)+Dtilde_c(i*2-1, g) + hRmv(i, g)*h_cmfd(i)
        do gp=1, ng
            if (gp .ne. g) mtfull(i+(gp-1)*totm_cmfd, l) = -hmSct(g, gp, i)*h_cmfd(i)
        end do
    end do
    mtflu = mtfull
    call LU_factorize(totm_cmfd*NG, mtflu)
end subroutine

subroutine matrix_shift(ke)
    use LU_general
    implicit none    
    
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "hmgz.FI"

    real(8), intent(in) :: ke
    integer :: i, l, m, g, gp
    real(8) :: ker
    ker = 1._8/ke
    mtflu = mtfull
    do g=1, ng
        do i=1, totm_cmfd
            l = (g-1)*totm_cmfd + i
            do gp=1, ng
                m = (gp-1)*totm_cmfd + i
                mtflu(m,l) = mtflu(m,l)-ker*Xchi(1,g)*hNuF(i,gp)*h_cmfd(i)
            end do
        end do    
    end do
    call LU_factorize(totm_cmfd*NG, mtflu)
end subroutine

subroutine makeQ_wldt(k, q)
    implicit none
    include "pspec.FI"
    include "arrays.FI"
    include "xsec.FI"
    real(8), intent(in) :: k
    real(8), dimension(:), intent(out) :: q
    integer :: i,  g, l, gp
    q =0
    do g=1, ng
        do i=1, totm_cmfd
            l = (g-1)*totm_cmfd + i
            q(l) = Xchi(1,g)*psi(i)/k
        end do
    end do
end subroutine


subroutine fdm_coarsemesh(initk)
    use define
    use lu_tdiag
    use powermethod
    implicit none
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "files.FI"
    real(8), intent(out) :: initk
    real(8) :: k, err, dr
    logical :: op
    integer :: iter
    logical :: add, wnd
    op = FALSE
    add = .false.
    wnd = .false.
!    call condense(add, wnd, .false.)
    call evaluateDtilde
    call make_matrix_withoutcurrent
    k=1.
    call powermethodG(totm_cmfd, 3, mt_cmfd, phi, psi, makeFS, makeQ, solve_LU, ccpsi, k, err, dr, iter)
    
    psi=psi*totm_cmfd/sum(psi)
    
    write(*, '(2(f12.5))'), k, dr
    write(hfiter, '(2(f12.5))'), k, dr
    initk = k
!    call testc_1g(k)
end subroutine

subroutine cmfd_EDT(k, err, iter, lmwnd)
    use define
    use lu_general
    use powermethod
    implicit none
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "files.FI"
    include "hmgz.FI"

    real(8) :: k, err, rerr(3), dr
    logical :: op
    integer :: iter, i, g
    real(8) :: ratio, rnorm
    logical :: lmwnd, add
    real(8) :: relfl1(totm_cmfd), relfl2(totm_cmfd-1)

    op = FALSE

!    phi_f = 1.
    if (NMVWND > 0) then
        add = .true.
    else
        add = .false.
    endif
!    call condense(add, lmwnd, .true.)
    call evaluateDtilde

    !phi_cp = phi_c
    call evaluateDhat
    !call make_matrix
    call make_matrix_full

!    ratio =sum(phi_cp)/sum(phi_c)
!    phi_c = phi_c*ratio
    
!    call powermethod_wldt(totm_cmfd, mtflu, makeFS, makeQ_wldt, matrix_shift, ccpsi, phi_c, psi_c, k, err, dr, iter)
    call powermethodW(totm_cmfd, totm_cmfd, mtflu, phi, psi, makeFS, makeQ_wldt, matrix_shift, solve_LU, ccpsi, k, err, dr, iter)
!    call prolongation
    
    if (PDFREF) then
        call err_coarse(phi, rerr(1), rerr(2), rerr(3))
        write(*, '(2(2x,i5,2x,f7.5,2x,e10.3,e10.3,e10.3,e10.3))',advance="no"), iter, k, err, rerr(1),rerr(2),rerr(3)
        write(hfiter, '(2(2x,i5,2x,f7.5,2x,e10.3,e10.3,e10.3,e10.3))',advance="no"), iter, k, err, rerr(1),rerr(2),rerr(3)
    else
        write(*, '(2(2x,i5,2x,f7.5,2x,e9.3))',advance="no"), iter, k, err
        write(hfiter, '(2(2x,i5,2x,f7.5,2x,e9.3))',advance="no"), iter, k, err
    endif    
    call makeFS
    
    rnorm = 1._8/sum(psi)
    phi=phi*rnorm
    write(fcmfd_EDT, '(2000(e11.5,1x))') k,(phi(i,1),i=1,totm_cmfd)
#ifdef _COEFF
!    write(fcp, '(2000(e11.5,1x))') ((hphi(i,g),i=1,totm_cmfd), g=1, ng)
!    write(fcj, '(2000(e11.5,1x))') ((hj(i*2,g),i=1,totm_cmfd), g=1, ng)
!    write(fcj, '(2000(e11.5,1x))') ((mtfull(i,i),i=1,totm_cmfd), g=1, ng)
!    write(fMTCMFD,      '(2000(e11.5,1x))') (mtfull(i,i),i=1,totm_cmfd)
!    write(fMTCMFDOFF,   '(2000(e11.5,1x))') (mtfull(i-1,i),i=2,totm_cmfd)    
    i=1
    relfl1(i) =  -mtfull(i+1,i)/mtfull(i,i)
    do i=2, totm_cmfd-1
        relfl1(i) =-( mtfull(i-1,i) + mtfull(i+1,i) )/mtfull(i,i)
        relfl2(i) = -mtfull(i+1,i)/mtfull(i,i)
    end do
    relfl1(i) = -mtfull(i-1,i)/mtfull(i,i)
!    write(fRL1CMFD,   '(2000(e11.5,1x))') (relfl1(i),i=1,totm_cmfd)    
!    write(fRL2CMFD,   '(2000(e11.5,1x))') (relfl2(i),i=2,totm_cmfd-1)    
        
#endif    
end subroutine

subroutine testc_1g(k)
    implicit none
    include "pspec.FI"
    include "xsec.FI"
    include "arrays.FI"
    include "hmgz.FI"
    real(8), intent(in) :: k
    integer :: i,j, o, m
    real(8) :: leak(3), absr(3), src(3), val(3), tot(3)
    real(8) :: diff(300), diff2(150), phis, phim, phio
    ! 1 group test
    call makeCurrent
    do i=1, totm_cmfd
        leak(1) = hj(i*2-1,1)+hj(i*2,1)
        leak(2) = hj(i*2-1,2)+hj(i*2,2)
        leak(3) = leak(1) + leak(2)
        absr(1) = phi(i,1)*hRmv(i, 1)*h_cmfd(i)
        absr(2) = phi(i,2)*hRmv(i, 2)*h_cmfd(i)
        src (1) = 1._8/k*psi_c(i)
        src (2) = hmSct(2,1,i)*phi(i,1)*h_cmfd(i)
        src (3) = src(1)
        absr(3) = absr(1)+absr(2)-src(2)
        tot (1) = hTot(i,1)*phi(i,1)*h_cmfd(i)
        tot (2) = hTot(i,2)*phi(i,2)*h_cmfd(i)
        tot (3) = (tot(1)+tot(2))/(phi_f(i,1)+phi(i,2))/h_cmfd(i)
        val (1) = (leak(1) + absr(1) - src(1))/phi(i,1)
        val (2) = (leak(2) + absr(2) - src(2))/phi(i,2)
        val (3) = (leak(3) + absr(3) - src(3))/(phi(i,1)+phi(i,2))
        diff2(i) = 1._8/3/tot(3)
!        if (i<150) diff(i) = -(jf(i*2,1)+jf(i*2,2))/(phi_f(i+1,1)+phi_f(i+1,2)-(phi_f(i,1)+phi_f(i,2)))*h(i)
        j = i*2-1
        if (j>1 .and. j<2*totm_cmfd) then 
            m = i
            o = i-1
            phio = phi(o,1)+phi(o,2)
            phim = phi(m,1)+phi(m,2)
            phis = (sum(hDif(o,:))*phio + sum(hdif(m,:))*phim)/(sum(hDif(o,:))+sum(hDif(m,:)))
            diff(j) = -h_cmfd(i)/2*(jf(j,1) + jf(j,2))/(phis-phim)
        endif
    end do
    contains
    
subroutine makeCurrent
    use define
    implicit none
    include "pspec.FI"
    include "arrays.FI"
    include "xsec.FI"
    integer :: g, i, i_f, i_c
    
    ! surface current condense    
    do g=1, ng
        do i=1, totm_cmfd
            if (i.eq.1) then
!            hj(i*2-1, g) = Dtilde_c(i*2-1, g) * hphi(i, g)
            hj(i*2-1, g) = Dtilde_c(1, g) * hphi(i, g)
            hj(i*2, g) = -Dtilde_c(i*2, g) * (hphi(i+1, g) - hphi(i, g))
            elseif(i.eq.totm_cmfd) then
            hj(i*2-1, g) = -Dtilde_c(i*2-1, g) * (hphi(i-1, g) - hphi(i, g))
!            hj(i*2, g) = Dtilde_c(i*2, g) * hphi(i, g)
            hj(i*2, g) = Dtilde_c(1, g) * hphi(i, g)
            else
            hj(i*2-1, g) = -Dtilde_c(i*2-1, g) * (hphi(i-1, g) - hphi(i, g))
            hj(i*2, g) = -Dtilde_c(i*2, g) * (hphi(i+1, g) - hphi(i, g))
            end if 
        end do
    end do
end subroutine
    
end subroutine 
end module
